Ultrafest

1 Data

num.splines <- 42
first.columns <- c(
    "speaker",
    "seconds",
    "rec.date",
    "prompt",
    "label",
    "TT.displacement.sm",
    "TT.velocity",
    "TT.velocity.abs",
    "TD.displacement.sm",
    "TD.velocity",
    "TD.velocity.abs"
    )
columns <- c(first.columns,
             paste0(rep(c("X", "Y"), num.splines),
                    "_",
                    rep(1:num.splines, each = 2)
                    )
             )

splines.raw <- list.files(path = "./pilot/results",
                   pattern = "*-aaa.csv",
                   full.names = TRUE) %>%
    map_df(~read.delim(., col.names = columns, na.strings = "*"))

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

splines.closure.raw <- list.files(path = "./pilot/results",
                   pattern = "*-closure.csv",
                   full.names = TRUE) %>%
    map_df(~read.delim(., col.names = columns, na.strings = "*"))

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

splines.phantom.raw <- list.files(path = "./pilot/results",
                   pattern = "*-phantom.csv",
                   full.names = TRUE) %>%
    map_df(~read_tsv(., col_names = columns, na = "*"))

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )

## See spec(...) for full column specifications.

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )

## See spec(...) for full column specifications.

rm(num.splines, first.columns, columns)

languages <- read_csv("./pilot/stimuli/languages.csv")

## Parsed with column specification:
## cols(
##   speaker = col_character(),
##   language = col_character()
## )

words <- read_csv("./pilot/stimuli/nonce.csv")

## Parsed with column specification:
## cols(
##   item = col_integer(),
##   word = col_character(),
##   ipa = col_character(),
##   c1 = col_character(),
##   c1phonation = col_character(),
##   vowel = col_character(),
##   anteropost = col_character(),
##   height = col_character(),
##   c2 = col_character(),
##   c2phonation = col_character(),
##   c2place = col_character(),
##   language = col_character()
## )

splines <- splines.raw %>%
    gather(spline, coordinate, matches("[XY]_")) %>%
    separate(spline, c("axis", "fan"), convert = TRUE) %>%
    spread(axis, coordinate) %>%
    mutate(word = word(prompt, 2)) %>%
    left_join(y = languages) %>%
    left_join(y = words) %>%
    mutate_if(is.character, as.factor)

## Joining, by = "speaker"

## Joining, by = c("word", "language")

splines.it <- filter(splines, language == "italian")
splines.pl <- filter(splines, language == "polish")

closure <- splines.closure.raw %>%
    gather(spline, coordinate, matches("[XY]_")) %>%
    separate(spline, c("axis", "fan"), convert = TRUE) %>%
    spread(axis, coordinate) %>%
    mutate(word = word(prompt, 2)) %>%
    left_join(y = languages) %>%
    left_join(y = words) %>%
    mutate_if(is.character, as.factor)

## Joining, by = "speaker"
## Joining, by = c("word", "language")

closure.it <- filter(closure, language == "italian")
closure.pl <- filter(closure, language == "polish")

2 Italian

2.1 Plot at maximum constriction

filter(splines.it, label == "max_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 754 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "max_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 1127 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "max_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 754 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "max_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 1127 rows containing non-finite values (stat_smooth).

2.2 Plot at peak 1

filter(splines.it, label == "peak1_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 844 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "peak1_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 997 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "peak1_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 844 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "peak1_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 997 rows containing non-finite values (stat_smooth).

2.3 Plot at closure

filter(closure.it, c2place == "coronal") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 937 rows containing non-finite values (stat_smooth).

filter(closure.it, c2place == "velar") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 1150 rows containing non-finite values (stat_smooth).

filter(closure.it, c2place == "coronal") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 937 rows containing non-finite values (stat_smooth).

filter(closure.it, c2place == "velar") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 1150 rows containing non-finite values (stat_smooth).

2.4 GAMMs at maximum constriction

splines.it$c2phonation.ord <- as.ordered(splines.it$c2phonation)
splines.it$c2place.ord <- as.ordered(splines.it$c2place)
splines.it$vowel.ord <- as.ordered(splines.it$vowel)
contrasts(splines.it$c2phonation.ord) <- "contr.treatment"
contrasts(splines.it$c2place.ord) <- "contr.treatment"
contrasts(splines.it$vowel.ord) <- "contr.treatment"
max.it <- filter(splines.it, label %in% c("max_TT", "max_TD")) %>%
    na.omit()

max.it.gam <- bam(
    Y ~
        c2phonation.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr"),
    data = max.it,
    method = "ML"
)

summary(max.it.gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -4.8549     0.1839 -26.406  < 2e-16 ***
## c2phonation.ordvoiceless   1.3251     0.2592   5.113 3.27e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F p-value    
## s(X)                          8.271  8.772 774.019 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 1.944  2.464   4.707 0.00501 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.676   Deviance explained = 67.6%
## -ML =  21185  Scale est. = 95.775    n = 5721

max.it.gam.null <- bam(
    Y ~
        s(X, bs = "cr"),
    data = max.it,
    method = "ML"
)

compareML(max.it.gam, max.it.gam.null)

## max.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## max.it.gam.null: Y ~ s(X, bs = "cr")
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf  Chisq    Df   p.value Sig.
## 1 max.it.gam.null 21203.08   3                            
## 2      max.it.gam 21185.13   6 17.943 3.000 7.912e-08  ***
## 
## AIC difference: -31.72, model max.it.gam has lower AIC.

plot_smooth(max.it.gam, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * X : numeric predictor; with 30 values ranging from -44.040500 to 61.699900.

plot_diff(max.it.gam, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * X : numeric predictor; with 100 values ranging from -44.040500 to 61.699900.

## 
## X window(s) of significant difference(s):
##  -44.040500 - -2.385191

acf_plot(resid(max.it.gam), split_by=list(max.it$rec.date))

max.it.gam.2 <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.it,
    method = "fREML"
)

summary(max.it.gam.2)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -6.3151     0.2486 -25.399  < 2e-16 ***
## c2phonation.ordvoiceless   1.3276     0.2281   5.820 6.21e-09 ***
## c2place.ordvelar           2.3702     0.2340  10.128  < 2e-16 ***
## vowel.ordo                 1.8548     0.2777   6.679 2.64e-11 ***
## vowel.ordu                -0.2606     0.3050  -0.854    0.393    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F  p-value    
## s(X)                          7.728  8.371 414.093  < 2e-16 ***
## s(X):c2phonation.ordvoiceless 2.765  3.469   8.285 7.17e-06 ***
## s(X):c2place.ordvelar         8.112  8.642 102.386  < 2e-16 ***
## s(X):vowel.ordo               5.395  6.592  15.107  < 2e-16 ***
## s(X):vowel.ordu               8.666  8.956  30.062  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.75   Deviance explained = 75.2%
## fREML =  20473  Scale est. = 73.809    n = 5721

max.it.gam.null <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr"),
    data = max.it,
    method = "fREML"
)

compareML(max.it.gam.2, max.it.gam.null)

## max.it.gam.2: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## max.it.gam.null: Y ~ c2phonation.ord + c2place.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr") + s(X, by = c2place.ord, bs = "cr")
## 
## Chi-square test of fREML scores
## -----
##             Model    Score Edf   Chisq    Df  p.value Sig.
## 1 max.it.gam.null 20751.11   9                            
## 2    max.it.gam.2 20473.03  15 278.072 6.000  < 2e-16  ***
## 
## AIC difference: -588.81, model max.it.gam.2 has lower AIC.

plot_smooth(max.it.gam.2, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -44.040500 to 61.699900.

plot_diff(max.it.gam.2, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -44.040500 to 61.699900.

## 
## X window(s) of significant difference(s):
##  -44.040500 - -8.793700

max.it.gamm <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.it,
    method = "fREML"
)

## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.

summary(max.it.gamm)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, 
##     by = vowel.ord)
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -6.3692     1.6283  -3.911 9.29e-05 ***
## c2phonation.ordvoiceless   1.4895     0.4027   3.699 0.000218 ***
## c2place.ordvelar           2.3929     0.4242   5.641 1.77e-08 ***
## vowel.ordo                 1.7681     0.5068   3.489 0.000489 ***
## vowel.ordu                -0.2405     0.5289  -0.455 0.649397    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df      F  p-value    
## s(X)                            7.338   8.048 20.434  < 2e-16 ***
## s(X):c2phonation.ordvoiceless   2.520   3.103  2.951   0.0297 *  
## s(X,rec.date)                 261.634 895.000  0.783  < 2e-16 ***
## s(X,speaker)                    4.209   8.000 12.604  < 2e-16 ***
## s(X):c2place.ordvelar           8.243   8.709 86.897  < 2e-16 ***
## s(X):vowel.ordo                 5.596   6.768  7.128 2.89e-08 ***
## s(X):vowel.ordu                 8.675   8.949 23.533  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.784   Deviance explained = 79.6%
## fREML =  20286  Scale est. = 63.649    n = 5721

max.it.gam.null <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.it,
    method = "fREML"
)

## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.

compareML(max.it.gamm, max.it.gam.null)

## max.it.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, 
##     by = vowel.ord)
## 
## max.it.gam.null: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord)
## 
## Chi-square test of fREML scores
## -----
##             Model    Score Edf  Chisq    Df   p.value Sig.
## 1 max.it.gam.null 20320.94  18                            
## 2     max.it.gamm 20286.36  21 34.583 3.000 6.441e-15  ***
## 
## AIC difference: -69.31, model max.it.gamm has lower AIC.

plot_smooth(max.it.gamm, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -44.040500 to 61.699900. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:51:50. 
##  * speaker : factor; set to the value(s): sdc02.

plot_diff(max.it.gamm, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -44.040500 to 61.699900. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:51:50. 
##  * speaker : factor; set to the value(s): sdc02.

## 
## X window(s) of significant difference(s):
##  -44.040500 - -9.861785

acf_plot(resid(max.it.gamm), split_by=list(max.it$rec.date))

2.4.1 Autoregressive correction

max.it <- arrange(max.it, rec.date, fan) %>%
    mutate(rec.date = as.character(rec.date))
previous <- ""
for (i in 1:nrow(max.it)) {
    if (max.it$rec.date[i] != previous) {
        max.it$start.event[i] <- TRUE
        previous <- max.it$rec.date[i]
    } else {
        max.it$start.event[i] <- FALSE
        previous <- max.it$rec.date[i]
    }
}

r1 <- start_value_rho(max.it.gam.2)

max.it.gam.AR <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.it,
    method = "fREML",
    rho = r1,
    AR.start = max.it$start.event
)

summary(max.it.gam.AR)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -6.3195     0.2587 -24.428  < 2e-16 ***
## c2phonation.ordvoiceless   1.3261     0.2367   5.602 2.22e-08 ***
## c2place.ordvelar           2.3922     0.2416   9.902  < 2e-16 ***
## vowel.ordo                 1.7830     0.2891   6.168 7.39e-10 ***
## vowel.ordu                -0.5921     0.3083  -1.920   0.0549 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F  p-value    
## s(X)                          7.950  8.534 398.352  < 2e-16 ***
## s(X):c2phonation.ordvoiceless 2.811  3.525   7.105 4.24e-05 ***
## s(X):c2place.ordvelar         8.085  8.631  96.067  < 2e-16 ***
## s(X):vowel.ordo               5.491  6.715  13.787  < 2e-16 ***
## s(X):vowel.ordu               8.673  8.958  27.925  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.749   Deviance explained = 75.1%
## fREML =  18712  Scale est. = 43.739    n = 5721

acf_plot(resid(max.it.gam.AR), split_by=list(max.it$rec.date))

acf_resid(max.it.gam.AR, split_pred = "AR.start")

plot_smooth(max.it.gam.AR, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -44.040500 to 61.699900.

plot_diff(max.it.gam.AR, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -44.040500 to 61.699900.

## 
## X window(s) of significant difference(s):
##  -44.040500 - -7.725615

2.5 GAMMs at closure

closure.it$c2phonation.ord <- as.ordered(closure.it$c2phonation)
closure.it$c2place.ord <- as.ordered(closure.it$c2place)
closure.it$vowel.ord <- as.ordered(closure.it$vowel)
contrasts(closure.it$c2phonation.ord) <- "contr.treatment"
contrasts(closure.it$c2place.ord) <- "contr.treatment"
contrasts(closure.it$vowel.ord) <- "contr.treatment"
closure.it <- closure.it %>%
    na.omit()

closure.it.gam <- bam(
    Y ~
        c2phonation.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr"),
    data = closure.it,
    method = "ML"
)

summary(closure.it.gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                2.8198     0.1835  15.371  < 2e-16 ***
## c2phonation.ordvoiceless   1.0034     0.2616   3.836 0.000126 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df      F  p-value    
## s(X)                          8.158  8.725 584.97  < 2e-16 ***
## s(X):c2phonation.ordvoiceless 2.819  3.514   7.04 4.17e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.514   Deviance explained = 51.4%
## -ML =  33952  Scale est. = 147.58    n = 8665

closure.it.gam.null <- bam(
    Y ~
        s(X, bs = "cr"),
    data = closure.it,
    method = "ML"
)

compareML(closure.it.gam, closure.it.gam.null)

## closure.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## closure.it.gam.null: Y ~ s(X, bs = "cr")
## 
## Chi-square test of ML scores
## -----
##                 Model    Score Edf  Chisq    Df   p.value Sig.
## 1 closure.it.gam.null 33969.15   3                            
## 2      closure.it.gam 33951.84   6 17.307 3.000 1.470e-07  ***
## 
## AIC difference: -33.44, model closure.it.gam has lower AIC.

plot_smooth(closure.it.gam, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * X : numeric predictor; with 30 values ranging from -47.595900 to 64.280100.

plot_diff(closure.it.gam, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * X : numeric predictor; with 100 values ranging from -47.595900 to 64.280100.

## 
## X window(s) of significant difference(s):
##  -47.595900 - -15.954203

acf_plot(resid(closure.it.gam), split_by=list(closure.it$rec.date))

closure.it.gam.2 <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = closure.it,
    method = "fREML"
)

summary(closure.it.gam.2)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                2.3530     0.2688   8.754  < 2e-16 ***
## c2phonation.ordvoiceless   1.2010     0.2437   4.929 8.43e-07 ***
## c2place.ordvelar           0.7587     0.2475   3.066  0.00218 ** 
## vowel.ordo                 1.9912     0.3036   6.559 5.72e-11 ***
## vowel.ordu                -1.9065     0.3386  -5.630 1.86e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df      F  p-value    
## s(X)                          6.801  7.659 329.07  < 2e-16 ***
## s(X):c2phonation.ordvoiceless 3.186  3.941  12.08 1.33e-09 ***
## s(X):c2place.ordvelar         7.512  8.274  75.43  < 2e-16 ***
## s(X):vowel.ordo               6.598  7.719  20.34  < 2e-16 ***
## s(X):vowel.ordu               8.376  8.862  25.00  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.58   Deviance explained = 58.2%
## fREML =  33340  Scale est. = 127.38    n = 8665

closure.it.gam.null <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr"),
    data = closure.it,
    method = "fREML"
)

compareML(closure.it.gam.2, closure.it.gam.null)

## closure.it.gam.2: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## closure.it.gam.null: Y ~ c2phonation.ord + c2place.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr") + s(X, by = c2place.ord, bs = "cr")
## 
## Chi-square test of fREML scores
## -----
##                 Model    Score Edf   Chisq    Df  p.value Sig.
## 1 closure.it.gam.null 33675.81   9                            
## 2    closure.it.gam.2 33339.51  15 336.302 6.000  < 2e-16  ***
## 
## AIC difference: -692.69, model closure.it.gam.2 has lower AIC.

plot_smooth(closure.it.gam.2, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -47.595900 to 64.280100.

plot_diff(closure.it.gam.2, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -47.595900 to 64.280100.

## 
## X window(s) of significant difference(s):
##  -47.595900 - -13.694082

closure.it.gamm <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
#        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = closure.it,
    method = "fREML"
)

## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.

summary(closure.it.gamm)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord)
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                1.8955     1.0209   1.857   0.0634 .
## c2phonation.ordvoiceless   1.4241     0.9217   1.545   0.1224  
## c2place.ordvelar           1.1709     0.9258   1.265   0.2060  
## vowel.ordo                 2.4567     1.1319   2.170   0.0300 *
## vowel.ordu                -1.8838     1.1323  -1.664   0.0962 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf   Ref.df       F  p-value    
## s(X)                            7.340    8.094 216.753  < 2e-16 ***
## s(X):c2phonation.ordvoiceless   3.641    4.424   5.772 6.97e-05 ***
## s(X,rec.date)                 535.105 1255.000   7.295  < 2e-16 ***
## s(X):c2place.ordvelar           8.082    8.631  83.016  < 2e-16 ***
## s(X):vowel.ordo                 7.253    8.233  13.469  < 2e-16 ***
## s(X):vowel.ordu                 8.327    8.826  32.646  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.796   Deviance explained =   81%
## fREML =  30882  Scale est. = 61.807    n = 8665

closure.it.gam.null <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
#        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = closure.it,
    method = "fREML"
)

compareML(closure.it.gamm, closure.it.gam.null)

## closure.it.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord)
## 
## closure.it.gam.null: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Chi-square test of fREML scores
## -----
##                 Model    Score Edf    Chisq    Df  p.value Sig.
## 1 closure.it.gam.null 33339.51  15                             
## 2     closure.it.gamm 30882.43  18 2457.080 3.000  < 2e-16  ***
## 
## AIC difference: -5746.74, model closure.it.gamm has lower AIC.

plot_smooth(closure.it.gamm, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -47.595900 to 64.280100. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:51:50.

plot_diff(closure.it.gamm, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -47.595900 to 64.280100. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:51:50.

## 
## X window(s) of significant difference(s):
##  -47.595900 - -19.344385

acf_plot(resid(max.it.gamm), split_by=list(max.it$rec.date))

2.5.1 Autoregressive correction

max.it <- arrange(max.it, rec.date, fan) %>%
    mutate(rec.date = as.character(rec.date))
previous <- ""
for (i in 1:nrow(max.it)) {
    if (max.it$rec.date[i] != previous) {
        max.it$start.event[i] <- TRUE
        previous <- max.it$rec.date[i]
    } else {
        max.it$start.event[i] <- FALSE
        previous <- max.it$rec.date[i]
    }
}

r1 <- start_value_rho(max.it.gam.2)

max.it.gam.AR <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.it,
    method = "fREML",
    rho = r1,
    AR.start = max.it$start.event
)

summary(max.it.gam.AR)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -6.3195     0.2587 -24.428  < 2e-16 ***
## c2phonation.ordvoiceless   1.3261     0.2367   5.602 2.22e-08 ***
## c2place.ordvelar           2.3922     0.2416   9.902  < 2e-16 ***
## vowel.ordo                 1.7830     0.2891   6.168 7.39e-10 ***
## vowel.ordu                -0.5921     0.3083  -1.920   0.0549 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F  p-value    
## s(X)                          7.950  8.534 398.352  < 2e-16 ***
## s(X):c2phonation.ordvoiceless 2.811  3.525   7.105 4.24e-05 ***
## s(X):c2place.ordvelar         8.085  8.631  96.067  < 2e-16 ***
## s(X):vowel.ordo               5.491  6.715  13.787  < 2e-16 ***
## s(X):vowel.ordu               8.673  8.958  27.925  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.749   Deviance explained = 75.1%
## fREML =  18712  Scale est. = 43.739    n = 5721

acf_plot(resid(max.it.gam.AR), split_by=list(max.it$rec.date))

acf_resid(max.it.gam.AR, split_pred = "AR.start")

plot_smooth(max.it.gam.AR, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -44.040500 to 61.699900.

plot_diff(max.it.gam.AR, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -44.040500 to 61.699900.

## 
## X window(s) of significant difference(s):
##  -44.040500 - -7.725615

2.6 GAMMs at peak 1

splines.it$c2phonation.ord <- as.ordered(splines.it$c2phonation)
contrasts(splines.it$c2phonation.ord) <- "contr.treatment"
peak.it <- filter(splines.it, label %in% c("peak1_TT", "peak1_TD"))

peak.it.gam <- bam(
    Y ~
        c2phonation.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr"),
    data = peak.it,
    method = "ML"
)

summary(peak.it.gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -5.8024     0.1809 -32.083   <2e-16 ***
## c2phonation.ordvoiceless   0.5019     0.2541   1.976   0.0482 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F p-value    
## s(X)                          8.691  8.958 562.990  <2e-16 ***
## s(X):c2phonation.ordvoiceless 1.695  2.137   0.932   0.425    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.59   Deviance explained =   59%
## -ML =  22127  Scale est. = 96.198    n = 5971

peak.it.gam.null <- bam(
    Y ~
        s(X, bs = "cr"),
    data = peak.it,
    method = "ML"
)

compareML(peak.it.gam, peak.it.gam.null)

## peak.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## peak.it.gam.null: Y ~ s(X, bs = "cr")
## 
## Chi-square test of ML scores
## -----
##              Model    Score Edf Chisq    Df p.value Sig.
## 1 peak.it.gam.null 22129.00   3                         
## 2      peak.it.gam 22126.83   6 2.171 3.000   0.227     
## 
## AIC difference: -0.63, model peak.it.gam has lower AIC.

## Warning in compareML(peak.it.gam, peak.it.gam.null): Only small difference in ML...

plot_smooth(peak.it.gam, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * X : numeric predictor; with 30 values ranging from -48.091200 to 64.315100.

plot_diff(peak.it.gam, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * X : numeric predictor; with 100 values ranging from -48.091200 to 64.315100.

## 
## Difference is not significant.

3 Polish

3.1 Plot at maximum constriction

filter(splines.pl, label == "max_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 568 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "max_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 635 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "max_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 568 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "max_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 635 rows containing non-finite values (stat_smooth).

3.2 Plot at peak 1

filter(splines.pl, label == "peak1_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 653 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "peak1_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 586 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "peak1_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 653 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "peak1_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)

## Warning: Removed 586 rows containing non-finite values (stat_smooth).

3.3 GAMMs at maximum constriction

splines.pl$c2phonation.ord <- as.ordered(splines.pl$c2phonation)
splines.pl$c2place.ord <- as.ordered(splines.pl$c2place)
splines.pl$vowel.ord <- as.ordered(splines.pl$vowel)
contrasts(splines.pl$c2phonation.ord) <- "contr.treatment"
contrasts(splines.pl$c2place.ord) <- "contr.treatment"
contrasts(splines.pl$vowel.ord) <- "contr.treatment"
max.pl <- filter(splines.pl, label %in% c("max_TT", "max_TD")) %>%
    na.omit()

max.pl.gam <- bam(
    Y ~
        c2phonation.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr"),
    data = max.pl,
    method = "ML"
)

summary(max.pl.gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              -8.47068    0.21100 -40.145   <2e-16 ***
## c2phonation.ordvoiceless  0.07766    0.30054   0.258    0.796    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F p-value    
## s(X)                          8.711  8.971 275.169  <2e-16 ***
## s(X):c2phonation.ordvoiceless 1.007  1.014   4.887  0.0271 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.347   Deviance explained = 34.8%
## -ML =  21329  Scale est. = 125.09    n = 5559

max.pl.gam.null <- bam(
    Y ~
        s(X, bs = "cr"),
    data = max.pl,
    method = "ML"
)

compareML(max.pl.gam, max.pl.gam.null)

## max.pl.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## max.pl.gam.null: Y ~ s(X, bs = "cr")
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df p.value Sig.
## 1 max.pl.gam.null 21331.36   3                         
## 2      max.pl.gam 21328.86   6 2.495 3.000   0.173     
## 
## AIC difference: -1.01, model max.pl.gam has lower AIC.

## Warning in compareML(max.pl.gam, max.pl.gam.null): Only small difference in ML...

plot_smooth(max.pl.gam, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * X : numeric predictor; with 30 values ranging from -36.642400 to 70.472300.

plot_diff(max.pl.gam, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * X : numeric predictor; with 100 values ranging from -36.642400 to 70.472300.

## 
## X window(s) of significant difference(s):
##  42.341167 - 70.472300

max.pl.gam.2 <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.pl,
    method = "fREML"
)

summary(max.pl.gam.2)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -8.2717     0.3101 -26.673   <2e-16 ***
## c2phonation.ordvoiceless  -0.1298     0.2776  -0.467   0.6402    
## c2place.ordvelar          -0.3388     0.2823  -1.200   0.2301    
## vowel.ordo                 0.5683     0.3433   1.655   0.0979 .  
## vowel.ordu                 0.6603     0.3387   1.950   0.0513 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F  p-value    
## s(X)                          8.625  8.917 102.950  < 2e-16 ***
## s(X):c2phonation.ordvoiceless 1.000  1.000   8.750  0.00311 ** 
## s(X):c2place.ordvelar         7.998  8.675  64.664  < 2e-16 ***
## s(X):vowel.ordo               3.554  4.512   5.854 5.38e-05 ***
## s(X):vowel.ordu               8.843  8.988  36.809  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.445   Deviance explained = 44.8%
## fREML =  20908  Scale est. = 106.33    n = 5559

max.pl.gam.null <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr"),
    data = max.pl,
    method = "fREML"
)

compareML(max.pl.gam.2, max.pl.gam.null)

## max.pl.gam.2: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## max.pl.gam.null: Y ~ c2phonation.ord + c2place.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr") + s(X, by = c2place.ord, bs = "cr")
## 
## Chi-square test of fREML scores
## -----
##             Model    Score Edf   Chisq    Df  p.value Sig.
## 1 max.pl.gam.null 21071.83   9                            
## 2    max.pl.gam.2 20907.82  15 164.019 6.000  < 2e-16  ***
## 
## AIC difference: -363.57, model max.pl.gam.2 has lower AIC.

plot_smooth(max.pl.gam.2, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): u. 
##  * X : numeric predictor; with 30 values ranging from -36.642400 to 70.472300.

plot_diff(max.pl.gam.2, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): u. 
##  * X : numeric predictor; with 100 values ranging from -36.642400 to 70.472300.

## 
## X window(s) of significant difference(s):
##  -36.642400 - -2.019467
##  38.013300 - 70.472300

max.pl.gamm <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1) +
        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.pl,
    method = "fREML"
)

## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.

summary(max.pl.gamm)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1) + s(X, speaker, bs = "fs", xt = "cr", m = 1, 
##     k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord)
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)               -5.9695     8.7701  -0.681   0.4961  
## c2phonation.ordvoiceless   0.7123     0.3220   2.212   0.0270 *
## c2place.ordvelar          -0.3770     0.3778  -0.998   0.3184  
## vowel.ordo                 0.7774     0.4645   1.674   0.0942 .
## vowel.ordu                 0.6505     0.4582   1.420   0.1558  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf   Ref.df       F  p-value    
## s(X)                            8.710    8.867  21.319  < 2e-16 ***
## s(X):c2phonation.ordvoiceless   1.001    1.001   0.042    0.838    
## s(X,rec.date)                 410.036 1587.000   0.869  < 2e-16 ***
## s(X,speaker)                    4.768    8.000 341.145  < 2e-16 ***
## s(X):c2place.ordvelar           8.136    8.710  58.659  < 2e-16 ***
## s(X):vowel.ordo                 7.917    8.606   5.720 7.04e-08 ***
## s(X):vowel.ordu                 8.777    8.952  25.767  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.805   Deviance explained = 82.1%
## fREML =  18415  Scale est. = 37.262    n = 5559

max.pl.gam.null <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.pl,
    method = "fREML"
)

## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.

compareML(max.pl.gamm, max.pl.gam.null)

## max.pl.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1) + s(X, speaker, bs = "fs", xt = "cr", m = 1, 
##     k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord)
## 
## max.pl.gam.null: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord)
## 
## Chi-square test of fREML scores
## -----
##             Model    Score Edf   Chisq    Df  p.value Sig.
## 1 max.pl.gam.null 18736.91  18                            
## 2     max.pl.gamm 18415.12  21 321.788 3.000  < 2e-16  ***
## 
## AIC difference: -82.32, model max.pl.gamm has lower AIC.

plot_smooth(max.pl.gamm, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)

## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): u. 
##  * X : numeric predictor; with 30 values ranging from -36.642400 to 70.472300. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:29:58. 
##  * speaker : factor; set to the value(s): ps02.

plot_diff(max.pl.gamm, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))

## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): u. 
##  * X : numeric predictor; with 100 values ranging from -36.642400 to 70.472300. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:29:58. 
##  * speaker : factor; set to the value(s): ps02.

## 
## X window(s) of significant difference(s):
##  4.472333 - 25.029700